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ABSTRACT 



Context. In spite of all the advances in multi-dimensional hydrodynamics, investigations of stellar evolution and stellar pulsations still 
depend on one-dimensional computations. This paper devises an alternative to the mixing-length theory or turbulence models usually 
adopted in modelling convective transport in such studies. 

Aims. The present work attempts to develop a time-dependent description of convection, which reflects the essential physics of con- 
vection and that is only moderately dependent on numerical parameters and far less time consuming than existing multi-dimensional 
hydrodynamics computations. 

Methods. Assuming that the most extensive convective patterns generate the majority of convective transport, the convective velocity 
field is described using two parallel, radial columns to represent up- and downstream flows. Horizontal exchange, in the form of fluid 
flow and radiation, over their connecting interface couples the two columns and allows a simple circulating motion. The main param- 
eters of this convective description have straightforward geometrical meanings, namely the diameter of the columns (corresponding 
to the size of the convective cells) and the ratio of the cross-section between up- and downdrafts. For this geometrical setup, the 
time-dependent solution of the equations of radiation hydrodynamics is computed from an implicit scheme that has the advantage 
of being unaffected by the Courant-Friedrichs-Lewy time-step limit. This implementation is part of the TAPIR-Code (short for The 
adaptive, implicit RHD-Code). 

Results. To demonstrate the approach, results for convection zones in Cepheids are presented. The convective energy transport and 
convective velocities agree with expectations for Cepheids and the scheme reproduces both the kinetic energy flux and convective 
overshoot. A study of the parameter influence shows that the type of solution derived for these stars is in fact fairly robust with respect 
to the constitutive numerical parameters. 

Key words, hydrodynamics — Cepheids — convection — methods: numerical 



1. Introduction 

Convection is one of the persistent problems in stellar astro- 
physics. Almost all stars contain regions where convective trans- 
port is important; in the photosphere, in the envelope, or in the 
interior where nuclear burning occurs. A description of convec- 
tion is therefore an essential ingredient to all types of investi- 
gations of stellar structure and evolution. Unfortunately, due to 
its nonlinear, nonlocal, and multi-length-scale nature, modelling 
convection turns out to be an intricate problem. 

Conceptually, there are several different approaches to the 
numerical simulation of convective transport in stars. The most 
straightforward approach is the time-dependent solution of the 
equations of radiation hydrodynamics in a 3D (or at least 2D) do- 
main to compute the convective flow patterns directly. However, 
this process is very expensive in computing time, prohibitively 
so for some applications. This particularly applies to problems 
with a large difference in relevant timescales, for instance be- 
tween the thermal and acoustic timescale in the stellar interior, or 
the hydrodynamics and radiative timescale in the outer layers of 
luminous stars. An additional limit is imposed by the restricted 
spatial resolution. Even the most elaborate, high-resolution sim- 
ulations of stellar convection are only capable of resolving the 
largest scales in the convective velocity field. The effect of turbu- 



lence on smaller length scales is effectively ignored, even though 
it is a possible interpretation to attribute the intrinsic numerical 
dissipation of the scheme to some (unknown) unresolved turbu- 
lence. In particular, hydrodynamics codes often include artificial 
viscosity for numerical stability, and 'unresolved turbulence' is 
the only physical mechanism that could be used to justify the 
inclusion of this additional dissipation. 

Despite the poor description of turbulence in hydrody- 
namics computa tions, multi-dim ensional simulations of solar 
granulation (e.g.lStein & Nordlun d. 1998; Asplund et aI.L l2000t 
IWedemever et aL , 20041) achieve remarkable quantitative agree- 
ment with observations. 

A completely different and more subtle approach than trying 
to resolve turbulence on large numerical grids are convection 
models that use an equation, or a set of equations, to describe 
convective transport either with a heuristic parametrization or 
based on turbulence theory. 

The most widely used of these I D-descriptions is the 
well known m ixing-length theory (MLT) (iBohm- Vitensel 1 1 9581; 
ICox & Giulii ri968). which originates in ideas of lPrandtlli 19251) . 
The MLT has been remarkably successful in stellar astrophysics 
in application to stars from white dwarfs to super giants, prob- 
ably because of its simple yet flexible parametrization and its 
robust reference to the adiabatic temperature gradient. 
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Modern alternatives to the MLT are convection models 
in which the original single-eddy assumption of the MLT 
has been replaced with a full spectruni of turbulence (e. 
ICanuto & MazzitelliL [19911 ICanutol [19961: ICanuto et all \l99' 
by either assuming or computing a turbulent energy spectrum. 
For some types of stars, these models can also avoid the mixing- 
length scale as a free parameter Convection models of this type 
are included in many stellar evolution codes as an alternative to 
MLT. 

However, all of these models are, in a similar way to the 
original MLT, local theories that do not provide any informa- 
tion about overshoot. In stellar-evolution codes, this deficit is 
overcome by adding overshoot by means of a separate (typically 
diffusive) parametrization. 

Another type of one-equation models consists of a time- 
dependent equatio r i for the turbulent kinetic energy (e 



S 'll99a 



LStellingwerfl 119821: iKuhfui 119861: iGehmevr & Winkle 
using heuristic approximations for individual terms. Containing 
a diffusion term for the turbulent kinetic energy, they also al- 
low a simple form of non-locality. These convection models are 
mainly g eared towards computatio ns of stellar pulsations ( e.g. 
Bono & Stellingwerft [l994l: iFeuchtingeil Il999bl: [ Kollath et al.l 



2002 1) where one is interested in the time-dependence of the con- 



vective transport. 

Finally, the most complete way to model turbulence is the 
'Reynolds stress approach' by so lving a set of moment equa tions 
(ICanutoill99ill993lll997l:lxiong . 1989; Xiong et al.L[T997l) that 
are terminated by a closure model at the third or fourth order. 
These closures are based on either the quasi-normal approxima- 
tion for the fourth-order moments or use a parametrization with 
reference to measured data, hydrodynamics (large-eddy) sim- 
ulations, or concepts such as the 'plume model' (see below). 
Turbulence models of this type are able to describe convective 
transport in a time-dependent and non-local way. For applica- 
tions of the Re ynolds s tress m odel to stellar surface convec - 
tion zones, see jKuEkal ([T999I) : iKupka & Montgomervl (l2002h : 



[Montgomery & Kupkal (|2002|) 

The two-column scheme, presented in this paper, is in- 
between the above categories, combining a hydrodynamics sim- 
ulation of the convective fluid flow with a parametrized, prede- 
fined geometry of the flow patterns. This setup is almost as sim- 
ple as a ID description; it describes up- and downstream with 
two parallel radial columns, and fluid flow over an interface be- 
tween those two columns allows a basic circulating convective 
motion. The two-column model could therefore be regarded as 
a simplistic 2D hydrodynamics scheme, limiting the horizontal 
range of the grid to just two cells. The very coarse description 
of the convective velocity field effectively implies that the most 
extensive convective flow patterns generate the majority of the 
convective transport. However, this assumption does not differ 
significantly from what is assumed in multi-dimensional hydro- 
dynamics computations that are also unable to resolve the full 
spectrum of convective turbulence. Since multi-dimensional hy- 
drodynamics achieve, in spite of this limitation, a good agree- 
ment with observations, the scenario of macroscopic convective 
patterns with distinct up- and downstream regions and little sub- 
structure (as also observed in the solar granula tion) appears to 
be a sufficient description of the actual physics (iNordlund et al.L 
11997 ). The two-column approach may therefore also give rea- 
sonable results. 

In the two-column scheme, the convective flux is computed 
directly from hydrodynamics and not from a heuristic model. 
Although it is not without numerical parameters, there is no 
'mixing-length parameter', nor anything equivalent. The method 



is also intrinsically non-local and the thickness of convective re- 
gions and the amount of overshoot are obtained consistently. 

The basic idea of modelling convection using separate ra- 
dial stratifications for up- and downstream regions is in fact an 
old one. In the 1960's to 1970's, predating advances in com- 
puting power that enabled 2D and 3D hydrodynamics com- 
putations to become possible, several similar two- or multi- 
stream models were devised. From observations, the solar gran- 
ulation pattern appeared to be separable to almost distinct hot 
and cool areas; it was therefore a logical first step to place two 
stratifications next to each other in order t o construct more re- 
alistic models of the s olar photosphere (Ma rgrave & Swih^3, 
19691: iNordlundl 1 19761). More recen tly a two-stream model 
was applied by iLesaffre et al.l (l2005h to investigate the con- 
vective Urea process in supernova-progenitor white dwarfs. 
In geophysics, a sim ilar concept , known as 'plume model', 
was introduced by [Morton et al.l (Il956l) . In its basic form, 
this model considers only plumes that are immersed in a 
static, surrounding medium, although there are also models 
that consider both up- and downdrafts and their interaction 
(e^Telford^ 1970; Wang & Albrecht, 1986; Chatfield & Brosl 
[l987[: [Randall et al.l [ 1 992[) . The idea of separated up- and down- 
wards streams was also used to construct closures for turbu- 
lence models (Abdella & McFarlane, 1997; Zilitinkevic h et al.1, 
1999; Lappen & Randan, 2001; Grvanik & Hartmanni [20021; 
Gryanik et al., 2005; Canuto et al., 2007) 

However, all existing multi-stream models differ from the 
present attempt in that the 'two-column-scheme' introduced be- 
low is based on fully implicit time-dependent radiation hydro- 
dynamics in both radial and horizontal directions without any 
ad-hoc assumptions or parametrizations for the physical cou- 
pling of the two columns. The term 'two-column' (in contrast to 
'two-stream') was chosen intentionally because the discretiza- 
tion scheme resembles that of two ID discretizations placed be- 
side each other in two parallel columns. In analogy to '2D' for 
'two-dimensional', we will occasionally refer to 'two-column' 
as '2C' in the following sections. 

The remaining paper is structured as follows. The next sec- 
tion. Sect. [2] introduces the two-column discretization scheme 
and its geometric derivation. The equations of radiation hydro- 
dynamics are given in Sect. [3] in analytical form, while Sect. [4| 
describes their discretization and considers details such as ar- 
tificial viscosity and radiative transport. Section [5] presents the 
deployed solution algorithm, followed by a demonstrating ex- 
ample and some parameter studies in Sect. [6] Finally, Sect. [7] 
draws the conclusions and summarizes the paper A verification 
of the method by comparison with detailed 2D hydrodynamics 
computations as well as applications in time-dependent calcula- 
tions of Cepheids' pulsations will be given in the forthcoming 
part II paper of this series. 



2. The two-column discretization scheme 

Figure [1] shows the setup of the two-column discretization 
scheme and the localization of the primary variables (for a com- 
plete listing of the primary variables, see Table[T]). The columns 
do not correspond directly to an individual convective cell but 
should be considered as a representation of all up- and down- 
stream flows, respectively. Correspondingly, the interface be- 
tween the two columns represents the sum of all contact surfaces 
between up- and downdrafts. 

In the 2C-scheme, the horizontal components of fluid flow 
and radiation, which actually occur in both the 6- and ^-direction 
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Fig. 1. The two-column discretization scheme and the localiza- 
tion of the primary variables (see Table [T]i. D is the typical dis- 
tance/diameter of the columns. The area shaded in gray repre- 
sents the discretization volume S_voli used for the scalar vari- 
ables pi, ei, Ji, and their respective equations. Advection oc- 
curs, as indicated by arrows, over the radial interface as well 
as over the interface Aintf between the two columns. The vector 
quantities m,-, Hr, ug, and Hg are included in their appropriate 
staggered-mesh location. Note that, although not illustrated in 
the figure, this all occurs in spherical geometry. 



Fig. 2. Spherical interpretation of the two-column scheme. The 
illustration on the left hand side shows the principle of how A' 
cone-like cells are assumed to be distributed over the sphere. 
Each of these cells, as sketched on the right, contributes to the 
interface area Aintf. Since the description is symmetric for the 
two columns (except for the relative cross-sections), these cones 
can be considered to represent either up- or downdrafts. 




The geometrical configuration of the discretization scheme 
is specified by two parameters. The first parameter is the typical 
horizontal length scale D, which can be interpreted as the di- 
ameter of the convective cells or the typical horizontal distance 
between up- and downdrafts. In contrast to the sketch in Fig.[Tl 
the two columns, in general, do not have the same size. The sec- 
ond parameter c/i (cf for column fraction) specifies the fraction 
of the sphere that is associated with column 1, and c/2 corre- 
spondingly, with c/2 = 1 - c/i, is that allocated to column 2. 

These two parameters, D and c/i , with their straightforward 
geometrical meaning are the main free parameters of the con- 
vection model. In Sect. 12.41 we will introduce a third constitutive 
parameter correlated with horizontal advection. Other numerical 
parameters, such as solution accuracy, radial advection scheme, 
artificial viscosity, boundary conditions, and grid resolution have 
only a minor effect on the solution. 

An equivalent yet more concrete quantity than D is the num- 
ber of convective cells on a sphere A^. In principle, it is possible 
to assign different values of A' to individual shells, or to spec- 
ify an analytical relation that defines N, for instance as a func- 
tion of the radius. However, in the absence of a robust physi- 
cal indication for the behavior of N, the present implementation 
uses the same N for every shell independent of the radius. That 
way, the up- and downdrafts are assumed to retain their identity 
throughout the convective region, which appears reasonable for 
photospheric convection zones of moderate depth. The relative 
cross-sections c/1 and c/2 must remain constant in all cases 
since changing their values would tilt the interface between the 
columns from the radial (coordinate) direction. 

By definition, D and N are related by 



D = 2r 



(1) 



which is based on the assumption of circular convective cells as 
sketched in Fig. |2] Since the formalism of the two columns is 
symmetric apart from in the relative cross-sections, it makes no 
difference whether these circular cells are considered as up- or 
as downdrafts. The right hand side of Fig. |2] illustrates the com- 
putation of the interface area between the two columns. Using 
Eq.[T]and summing for A columns, we obtain 



lintf - 



(2) 



Despite the geometric motivation given in Fig. |2l this picture 
should not be interpreted literally. By combining both the 6- 
and (^-directions of spherical geometry to one generic horizon- 
tal variability, the direct correlation with the three-dimensional 
configuration disappears. It is therefore not sensible to interpret 
expressions from the two-column formalism using, for example, 
a specific slice through the setup shown in Fig.|2] 

However, the only points where the geometrical configura- 
tion enters the scheme are the definitions of the interface area 
in Eq. |2]and of the horizontal derivatives in Eqs. |3]&|4] Since 
these equations are coupled with each other by the requirement 
to ensure that GauB's theorem applies also in the discrete case, a 
different geometrical picture would change Eqs.|2]-|4]by only a 
constant factor 



of spherical geometry, are described each by just one 'horizon- 
tal' variable. For these variables, ug and Hg respectively (see 
Table [TJ, the subscript '0' does not refer to spherical geometry 
components but is used more generally to denote 'horizontal' 
variables. 



Horizontal derivatives 

Using the typical horizontal length scale D, we can approximate 
derivatives in the horizontal direction by 

- " 7; ^"(^^ = T- VaV2 Ag(X), (3) 
r 06 D 2r 
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where '0' is again used to denote the generic horizontal exten- 
sion of the 2C-scheme. 

For second-order derivatives with respect to 6, we adopt the 
estimate 



Column 1 Column 2 



I d^X 2 



N 

N_ 
[ 4^2 



(X2 - Xi ) for column 1 



(Xi - X2) for column 2 



, (4) 



which is based on the assumption of a basically periodic varia- 
tion of X in 0-direction (and consequently of ^ as well) due to 
the alternating succession of up- and downdraft columns. This 
already requires over-stretching of the geometrical picture but 
since these derivatives exist at a less important point (in the ^ 
term of the viscous forces, Eqs. |40] & l44l i. the approximative 
evaluation of |^ is acceptable. 



2.1. Scalar discretization 

The discretization of scalar physical variables and equations 
uses discretization volumes similar to that highlighted in gray 
in Fig. [1] The volume of the scalar cells (distinguished by the 
prefix 'S_' for scalar) is computed to be the appropriate fraction 
of the shell between the radii r, and r,+i 







V_vol, 

.____(_________ 




Ur,i H„i 

















(5) 



S.voh - cfif{rl^-r^) 
S.V0I2 = cf2'f{rl^-r^) . 

The advective fluxes (transported 'volume' during a time step) 
for the scalar discretization are indicated by arrows in Fig. [1] 
Radial advection consists of two contributions; one from the 
proper motion of the fluid, and one due to movement of the adap- 
tive grid 



Fig. 3. Discretization volume for vector variables and equations 
in radial direction. The radial boundaries of the cell are defined 
by the averaged radii ?, and . The arrows illustrate horizontal 
advection which is composed of two parts correlated with ug, 
and Mfli-i- 



r,. Using the definition of r, we can now compute the volumes 
(with the prefix 'V_' for vector) of these cells 



V-Vol, = c/, i^i{rl^-rl^) 
V_vol2 = cf2'fi{rl^-rl^) 



(9) 



For the advective flux in the radial direction, we interpolate the 
velocity assuming flux conservation, i.e. 



S.fluxi = c/i [Att rj uij 6t-f (r ' 
S.flux2 = c/2 [47T rf U2J (Jf - y (r/ 



new3 
new'i 



i 



)] 



(6) 



= \[r^Ui + r]^^UM) 



'•,+l/2Mi+I/2 2 ■ 

and therefore in analogy with Eq. |6] we obtain 



(10) 



where 6t is the time step during which the grid adaptivity al- 



ters the radius of the grid point ; from r" to r"™ . Note that V_fluxi — cf\ 

the individual radial velocities uij and M2,i were used in the two 

columns. 

The horizontal advective flux between the two columns is 

computed to be V_flux2 = c/2 



TT I / 2 2^ 

S_fluxfl = Aintf ue,i6t = - ylNjl [r^^y - r,. j ue,i6t . 



(7) 



An\{ijui^i + rf^^ Mi,,+i) St 

T yi ~ ''i ) 

An \ [r] U2J + rf_^^ U2j+i) St 

3" yi ~ ''• ) 



(11) 



For both the horizontal velocity ug and the horizontal advective 
flux S_fluxe, a positive sign corresponds to a fluid flow from col- 
umn 1 to column 2 by convention. The analogous convention is 
also used for the horizontal radiative flux Hg. 

2.2. Vector discretization - radial 

Vector-type variables and equations are discretized on a stag- 
ge red mesh where th e radial part closely resembles that given 
bv lDorfietal] (l2006h . To define discretization volumes for the 
radial vector components, we start by defining 'averaged' radii 
located in between the grid point positions and denoted by ? 



-3 3 ^ f 3 , 3 \ 



(8) 



As indicated by the arrows in Fig. [3] the horizontal flux between 
the two radial vector volumes V_voli and V_vol2 is composed of 
two parts correlated with Mg, and Me,_i 

V.fluxe = ^ VaV2 [ifi - rf) ugj + {rf - ^h) Me,/-i] St . (12) 



2.3. Vector discretization - horizontal 

The discretization of the horizontal components of vector vari- 
ables and equations (namely of ug and Hg) uses discretization 
volumes as illustrated in Fig.|4] The corresponding volumes and 
fluxes are labeled with the prefix 'H_' for horizontal. The dis- 
cretization cell with H_vol is centered on the interface between 
the columns and considers half the volume of the sphere 



Figure[3]shows two of these averaged radii, r, and r,_i , which es- 
tablish the discretization volumes centered around the grid point 



H-V0l=iy (r^.i-r?) 



(13) 
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The second half of the volume is assumed to mirror the physics 
of the first one. Even though the discretization volume H_vol 
represents only one half of the shell, it therefore describes the 
horizontal components of the entire sphere. 

From comparison with Eq. |5] we observe that H_vol - 
1/2 (S_voli + S_vol2); the flux in the radial direction H_flux is 
assembled in the same way from S_fluxi and S_flux2 (see Eq.|6]) 



Column 1 Column 2 



H.flux = i 



A-nrf (cfiuij + cf2U2j) St 



(14) 



Since the discretization scheme represents a large number of 
convective cells distributed over the sphere, the two columns can 
be considered as part of a sequence of alternating up- and down- 
drafts. Along this sequence, the direction of horizontal fluid flow 
and radiation switches its sign repeatedly. This implies that a 
right-hand orientated flow (as illustrated in the lower part of 
Fig. lU is confronted with an equal flow in the opposing direc- 
tion when reaching the (in this case) right hand cell boundary. 
To allow for this effect, a dissipation term was included in the 
equation of motion that could be interpreted as annihilation of 
the momentum of the two opposing flows 



^'anhl = - |S_fluxH| (c/ipi -I- cf2P2)ug 



(15) 



where (c/ipi + C/2P2) ug is the momentum in the cell and 
|S_fluxH| describes the volume fraction swept against the hori- 
zontal boundary. In the equation of internal energy, this dissi- 
pated energy enters as 



^anhl 

= |S_fluxH| (c/ipi +cf2P2)u, 



e ■ 



(16) 



The location at which this energy is deposited depends on the 
direction of horizontal flow. In the case sketched in Fig. H] the 
dissipated energy would be deposited into column 2. 

As we will make use of it in the set of discrete equations 
(TableO, we finally define pg, the averaged density appropriate 
for the 'horizontal' discretization volume 



Pe = (c/iPi + C/2P2) 



(17) 



Despite the similar notation, this horizontally averaged density 
should not be mistaken for the radially averaged densities pj and 
P2 that are introduced in Sect. 14.51 

2.4. Horizontal advection and energy conservation 
Horizontal advection 



m g a sec ( 
l ll977h . 



but 



order van Leer advection scheme (I van Leetl 1 197 
for horizontal advection, i.e. fluid flow from one column to the 
other, a higher order scheme is obviously not applicable. The 
most straightforward approach is donor cell advection, but, in 
general, we can allow some variation within the columns, as il- 
lustrated in Fig.|5]for the example of density. This produces the 
following advection scheme 



j (1-A)pi+Ap2 for Col. 1 ^ Col. 2 
^"\(1-A)p2 + Api for Col. 1 ^ Col. 2 



(18) 



where we adopt the convention of denoting advected quantities 
with an overhead tilde. From Eq. [18] one recovers both simple 
donor cell advection for A - and centering between pi and 
P2 for A = 1/2. However, it is advisable to retain a small value 
of A, i.e. resulting in an advection scheme similar to donor-cell. 







H.vol 




He, 






1 


r 









Fig. 4. Discretization volume for the horizontal components of 
vector variables and equations. The radial advective flux is com- 
posed of two parts that, except for a factor 1 /2, resemble the 
scalar radial fluxes S_fluxi and S_flux2. The gray arrows in the 
lower part illustrate the mirroring principle used when consider- 
ing the two columns as part of a sequence of up- and downdrafts. 
In that picture, the two outer cells with the left-hand arrows are 
actually one and the same. 



Column 1 



Column 2 




Fig. 5. Reconstruction scheme for horizontal advection. 
Allowing some variation in the variables (here p as an exam- 
ple) across the width of the columns decreases the contrast at 
the interface between up- and downdraft. For computation of the 
horizontally advected quantity p, this variation is assumed to be 
proportional to the contrast between the two columns Ap. The 
parameter A with < A < 1/2 allows for a continuous transition 
between donor cell and centered advection. 



since centered advection can produce unrealistic values for ad- 
vected radial momentum: as part of a large circulating fluid flow, 
the convective motions in the two columns correspond to each 
other. Centering the momentum between the up- and downdraft 
column therefore gives almost zero momentum. Consequently, 
the up- and downdraft flows lose hardly any momentum due to 
horizontal advection, even if there is a large horizontal exchange 
of mass and internal energy. This mechanism only affects the 
momentum as it is the sole advected quantity where the values 
for the two columns usually have opposing signs. 

Despite this potentially unphysical behavior for larger val- 
ues of A, the formalism in Eq. [18] provides an additional free 
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parameter for adjusting the convection zones obtained from the 
2C-scheme with reference to established results. 

Energy conservation 

The equations of radiation hydrodynamics (Eqs. |26l- l3TTi are 
discretized conservatively. Hence, the scheme conserves mass, 
momentum, internal energy, as well as the moments of radi- 
ation. Although analytically equivalent, this does not translate 
into conservation of the total energy in the discrete case. Usually, 
this is not crucial for ID computations. Moreover, the adaptive 
grid provides a fine grid resolution for all gradients and accord- 
ingly minimizes spatial discretization errors. In case of the hori- 
zontal components in the 2C-scheme, we now have to consider a 
very coarse spatial representation where, in particular, advection 
from one column to another requires some attention. 

In the present discretization, advection of total energy con- 
sists of three components: internal, radiative, and kinetic energy. 
The former two are treated accurately by the advection terms in 
the corresponding equations of internal energy and radiation en- 
ergy. In contrast, the advection of kinetic energy is modeled only 
indirectly by density and momentum transport. Analytically, ad- 
vection of kinetic energy, momentum, and density are related by 

^■{u{pu^)^u V iupu)-^jU^V (up) . (19) 

Integration over a cell volume provides the discrete (approxi- 
mate) equivalent 

^ |pM?Flux,- ^ u ^ pui Flux; - jU^ ^ pi Flux/ (20) 

(' (' / 

where the quantities with a tilde are advected over / cell bound- 
aries with the transported volumes Flux,. Using this formula, we 
can now compute the kinetic energy effectively transported by 
the advection of mass and momentum. 

In application to horizontal advection from one column to 
the other, we obtain for column 1 

jpu^V Jluxe - uipu V_fluxe - ^Mjp VJluxg (21) 
and column 2 

jpu^V Jluxe - U2pu V_fluxe - jU^WJluxg . (22) 

Comparing these two lines, it becomes apparent that the ad- 
vected kinetic energy ^pu^ is not identical in both columns: a 
certain flux of density p and momentum pii over the interface 
between the two columns induces a change in kinetic energy 
in column 1 as given by Eq. [21] while column 2 experiences 
a change according to Eq. |22l The advection process therefore 
creates an error in the kinetic energy balance, and consequently 
also in the conservation of total energy. 

To allow for this deficit in the total energy balance, we com- 
pute the difference and place it as a source term into the equation 
of internal energy 

^adv = (mi - U2)'pu V_fluxfl --(«!- "2) P S_fluxe . (23) 

In the simplest case, horizontal transport uses donor cell advec- 
tion or, more generally, a formalism as in Eq.[T8] Assuming that 
we compute p and pu analogously, i.e. with the same A, we can 
further simplify 

£adv = (Ml - uifp' 5|V_fluxe| (24) 



where p* is given by 

r(l-^)pi-^P2 for V_fluxe>0 
^ \(1-A)p2-Api for V_fluxe<0 ^ ' 

and p is the radially averaged density (see Sect. 14.51 ). For A = Q 
(donor cell), p* becomes the upstream value, in which case we 
could write p* = p. 

From Eq. [24l we have E^dv > 0, i.e. E^dv always acts as a 
source term for the internal energy. E^^^ therefore effectively de- 
scribes the dissipation of kinetic energy in the course of advec- 
tion from one column to another This dissipation increases with 
the radial velocity difference \ui - M2I between the two columns. 
In a convection zone, the two columns hold opposing up- and 
downdraft motions and |mi - M2I is quite large; |mi - M2I - '^Umnv 
In these cases, the dissipation term becomes indispensable to the 
total energy balance; in the examples presented in Sect. |6] it can 
account for more than 30% of the energy throughput (i.e. lumi- 
nosity). 

Depending on the direction of the horizontal flow, the dissi- 
pated energy is deposited in the receiving column. In doing so, 
the contribution from Eq.|24]must be divided radially to be con- 
sistent with the scalar discretization of the equation of internal 
energy. 

2.5. Radial distribution of grid points 

In the preceding paragraphs, we constructed the two columns 
discretization scheme with reference to a given radial distribu- 
tion of the grid points r,. We now have to adopt a method to 
determine these grid point positions. 

For obvious reasons, the convective fluid flow prohibits a 
Lagrangian grid customarily used in stellar models. A spatially 
fixed Eulerian grid is also poorly suited to our needs for two rea- 
sons. First, advection alters the stellar structure. Starting from 
an initial, purely radiative model, the star shrinks significantly 
with the onset of advective transport. Secondly, this scheme is 
intended to be used in computing stellar pulsations, i.e. to fol- 
low the convective circulating motion while the entire envelope 
moves in- and outward in the course of stellar pulsation. 

To m eet these requirement s, the code uses an adaptive grid 
equation (lDorfi&DrurvilI987h that redistributes the grid points 
continuously according to the evolving physical structures and 
therefore provides high resolution as needed, e.g. at photo- 
spheric gradients, while following radial movements of these 
features due to structural changes or stellar pulsation. This adap- 
tive grid equation is solved implicitly together with the physi- 
cal equations. Since the grid equation is an elliptic differential 
equation, this approach is only possible with an implicit solving 
method. 

In the application to the two-column scheme, the same grid 
point distribution is used for both columns. Otherwise horizontal 
advection would become far more complicated and - a serious 
issue for the implicit solving algorithm - non-local with respect 
to the grid index /. 

Variables from both columns are used as 'grid- weights', in 
particular the grid adapts according to gradients in density, inter- 
nal energy, and V^d in each column. The grid resolution is there- 
fore increased in both columns in identical ways, even though, 
in general, only the physical structure in one of them actually 
demands this high grid resolution. 
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3. Physical equations 

The physics within the two-column geometrical setup is 
comput ed from t he equations of radiation hydrodynamics 
(e.g. Mi halas & M ihalas. 1984). The radiation field is thereby 
described using the first three gray moments of the intensity 
J, H, and K, which correspond to the radiative energy den- 
sity, radiative flux, and radiative pressure, respectively. An 
Eddington factor /edd closes the moment equations. Neglecting 
scattering, the source function of radiation, S, is given by 
the Stefan-Boltzmann law S - cr/nT*, and /cr and Kp are the 
Rosseland and Planck mean opacities. The gas pressure P 
and gas temperature T are given by the equation of state. Self 
gravity is described by the gravitational potential 0, which 
is assumed to be spherically symmetric, i.e. we do not allow 
for the (negligible) gravitational interaction between up- and 
downstreams. G is Newton's gravitational constant. Artificial 
viscosity, discussed in detail in Sect. 14.11 enters in the form of 
the viscous pressure tensor Q. 

The system of analytical equations is given by: 

Equation of continuity 

|-p + V-(Hp) = (26) 
ot 

Equation of motion 

—(pu) + V ■ (upu) H- VP H- pV^ - —K^H -H V • Q = (27) 
ot c 

Equation of internal energy 

Q 

—(pe) + V ■(upe) + PV u- 4kkpp(J - 5) -h Q : Vm = (28) 
ot 

Poisson equation 

Acf) = 47rGp (29) 
Radiation energy equation 

—J + 'V-{uJ) + cVH+K:Vu+c Kpp{J - 5) = (30) 
ot 

Radiation flux equation 

Q 

—H + V-(uH) + cV-K + HVu+cKRpH^O (31) 
at 

Radiation equations for liigli optical deptlis 

The difference (J-S) gradually vanishes with increasing optical 
depth, i.e. towards the interior of a star. 

Therefore, the coupling term (J-S) between radiative en- 
ergy and gas energy becomes numerically unresolvable for high 
optical depths. To derive still the correct contribution from this 
coupling for the equation of internal energy, the corresponding 
term Attkp p (J-S) is expressed by the radiation en ergy equation 
and in serted into the equation of internal energy (i Feuchtingeii 
Il999ah . This coiTesponds to evaluating the sum 'Equation of en- 
ergy' + ^ 'Radiationenergy equation', where terms with (7-5) 
cancel out each other. The conversion factor — relates the zeroth 

c 

moment of the intensity J to the radiation energy density 
-ipe+^J) + V-[u(pe+i^J)] + 

+PV ■u + ^K-.Vu+A-kV ■ H + Q:Vu =0 . (32) 



Inwards of a predefined stellar depth, the equation of internal 
energy is substituted with this sum, i.e. Eq.[32]is solved instead 
ofEq.|28] 



4. Discrete set of equations 

After introducing the analytical form of the equations of radi- 
ation hydrodynamics, we now develop their discrete version. 
Table [T] summarizes the primary variables, the coiTesponding 
discrete equations, and the closures of the system. As illustrated 
in Fig. [1] r and m, as well as the 'horizontal' variables ug and 
Hg, are integral quantities for both columns. All other variables, 
p, e, u, J, and H exist in duplicates, assigned individually to the 
two columns. 

4. 1 . Artificial viscosity 

In the continuum description of fluids, shock fronts - and, in 
the present case, horizontal shear flows - may become indef- 
initely sharp. In hydrodynamics codes, the smallest physical 
length scale is given by the mesh size of the numerical grid; 
on this length scale, numerical dissipation intrinsic to the spa- 
tial discretization becomes effective. In the present implemen- 
tation, the adaptive grid continuously refines to resolve all gra- 
dients properly on the grid. This reduces the intrinsic numeri- 
cal dissipation and can thus lead to a runaway effect of succes- 
sively steepening gradients and subsequent grid refinement. It is 
therefore necessary to include an artificial viscosity as a measure 
of broadening naiTow physical features on a predefined length 
scale. Consequently, this also limits the maximum grid resolu- 
tion to which the adaptive grid will be refined to. 

In this way, artificial viscosity, by specifying the minimum 
length scale in the computation, plays a more important role than 
in usual Lagrangian or Eulerian hydrodynamics codes. 

Due to the small overall dissipation of the numerical scheme, 
it is also sometimes necessary to include some extra viscos- 
ity to limit amplitudes and velocities, e.g. of stellar pulsations. 
However, in the results presented in Sect. |6] the influence of 
the artificial viscosity always remains negligible and is apparent 
only in a minor smoothing of velocity spikes. 

For the artificial vis cosity, the geometry-ind e pende nt de- 
scription provided by iTscharnuter & Winkled d 19791) was 
adopted. In this description, modeled by analogy with the or- 
dinary (molecular) fluid viscosity, the viscous pressure tensor 
reads 

Q = -fiQ J[VH],y™ - 1 i V ■ hJ (33) 

where the viscosity coefficient fiQ contains parameters for 'lin- 
ear' (pseudo-molecular) viscosity and 'quadratic' viscosity 
?quad (where 'quadratic' refers to the quadratic dependency on 
the velocity field, which causes it to act in a way similar to a 
turbulent viscosity) 

jUQ = ^lin^visc P Cj + ^quad^visc P (-V ■ H, 0) . (34) 

The use of the maximum implies that expanding flows are un- 
affected by viscosity. The viscous length scale /vise is set to 
the characteristic extension of the problem (and of the numer- 
ical grid), e.g. the radius in spherical geometry, c, is the local 
speed of sound. For the symmetric velocity gradient, the notation 
[Vulsym was introduced. The symmetric description ensures that 
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Table 1. Set of primary variables and the corresponding equations; for the discrete equations see also Table|3] 



Variable 


Description 


Equation 


ri 


Radius 


Adaptive grid equation 




llllCglalCQ IllaSo 


I^UlShUll CqUaLiUll, i.e. laUlal illLCgldLiUll Ui lllctaa 


pihPii 


Density 


Equation of continuity (in each column) 


PihPii 


Averaged density 


Radial averaging of p: Eq.[60]& Eo.l6T] 


en, 62/ 


Specific internal energy 


Equation of energy (in each column) 




Radial velocity 


Equation of motion, radial component (in each column) 




Horizontal velocity 


Equation of motion, horizontal component 


Jii, Jli 


0* moment of radiation 


Radiation energy equation (in each column) 




1-^' moment of radiation, radial 


Radiation flux Eq., radial component (in each column) 




P' moment of radiation, horizontal 


Radiation flux Eq., horizontal component 



Closures: - tabulated equation of state (temperature, gas pressure), evaluated separately in each column: 
Ti = r(pi,ei), T2 = T(p2,e2), Pi = P(pi,e,), P2 = P{p2,e2) 

- tabulated opacities (Rosseland mean), evaluated separately in each column: 

Kl = K(pi,e[),K2 = K(p2,e2) 

- closure of radiation moments with an Eddington factor /ejd = K/J = 1/3 



rotation, which does not affect the physical structure, remains 
unaffected by viscosity 



(35) 



The contributions of artificial viscosity to the equations of mo- 
tion and internal energy follow directly from the viscous pres- 
sure tensor. The viscous force is computed to be the divergence 
of the viscous pressure 



/q = V ■ Q 



(36) 



The viscous energy dissipation is obtained by contraction of the 
viscous pressure tensor with the gradient of the velocity field. 
Since Q is symmetric, there is no difference between using the 
velocity gradient Vh or the symmetric velocity gradient [VH]syn, 

eQ = Q : Vh . (37) 

To apply this recipe in the present case, Eqs.[33]-[32]must be 
evaluated by assuming spherical geometry. Since the 2C-scheme 
describes all types of horizontal variability and dynamics with 
only one interface between the two radial columns, the 9- and (p- 
directions of spherical coordinates are not considered separately 
and we adopt the identities ug = and ^ = To allow for 
that, all derivatives in the (;A-direction are assumed to be taken on 
the great circle, i.e. for 6 = n/l. Also note that Q is symmetric by 
definition and we therefore finally have four independent entries 
for the viscous pressure tensor in the two-column geometry: Qrr, 
Qr0, Qee, and Qg^. 

In principle, viscosity couples fluid flows in different coor- 
dinate directions. In the 2C-scheme, due to the combined dis- 
cretization of 9- and ^-components, the corresponding terms 
in the spherical symmetric description become ambiguous 
in interpretation; the two-column representation of the three- 
dimensional flow is too simplistic to enable a proper modelling 
of this effect. The viscous interaction between the two directions 
of fluid flow was therefore neglected by assuming -Q for the 
viscosity in the 0-direction, and ug - Q m the radial direction. 
For the viscosity in the radial direction, we then obtain 



Qre = 



2 / du,- Ur 



1 1 dUr 



(38) 
(39) 



f 3 5/3x2 dQrg 



do 



2 I dUr Ur 



For the viscosity in the ^-direction, we arrive at 



\ Idug Ug 

Qre = -/^Go ^ 

2\ dr r 

_ Wdug 



/Q« = 2^A(.^e.«) + 2i^(4e,,) 



rde 



dug Ug 



2 II dug 



(40) 
(41) 

(42) 
(43) 

(44) 
(45) 



In Eq|44l an additional factor 2 was included for fqg because 
we are considering forces in both the 6- and (^-directions, even 
though they are combined in the discretization process. 

In the discrete case, derivatives with respect to radius trans- 
form into differences between radial indices (Ar), and derivatives 
in the ^-direction are discretized using Eqs.|3]&|4] 

The various terms of the artificial viscosity (radial - horizon- 
tal, shear - non-shear) include separate coefficients juq to allow 
for their individual adjustment. Table |2] presents the jjq coeffi- 
cients with the parameters on which they depend. The computa- 
tion of the /iQ coefficients is similar to that described by Eq. |34] 
except for details related to the staggered-mesh location of the 
involved variables; the turbulent ('quadratic') viscosity param- 
eter is only used for the radial, non-shear part (/kqi and yUgi). 
Where appropriate, the jjq's are evaluated separately in each 
column, although the viscosity parameters are the same in both 
columns. In total, there are 5 viscosity parameters, although un- 
til now only three (except for testing purposes) were actually 
used in the computations. The default values adopted in the ex- 
amples presented in Sect. |6] are ^lin = 10"^, ^quad - 10"^, and 

^fishear - lO""*. 
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Table 2. Compilation of viscosity coefficients and parameters. 



Coefficient 


Direction of action 


Parameters 




radial 






radial 




/jQshear 


radial 


^shear 




horizontal 




^QBshear 


horizontal 


9eshear 



The discretization of the viscous force and energy dissipation 
uses the same discretization volumes as the corresponding equa- 
tions, i.e. the equation of motion and equation of internal en- 
ergy. To emphasize the volume-integrated variables, the discrete 
forces and energies are denoted by capital letters; Fq - J fq dV 

and Eq- J eq dV. 

For the viscous force in the radial direction, we obtain for 
column 1 



J7 f^^.i -il^yUl 



+jUQshearT-T ("l - Ul) -V_VOl 



(46) 



and for column 2 



Fq2 = -c/2-A,.U,r — 



N 



1, 



-jUQshearTT ("l " "2) -V_VOl . 

2 



(47) 



Note that both shear forces are discretized with jV.vol instead 
of V_voli and V_vol2 to allow them to cancel out each other for 
the two columns. 

The corresponding viscous energy dissipation reads 



£qi 



2 / ArMi 



3 \ Arr 



yl S_voli 



A'Qshear ^ 



N lu\ — U2^^ 



iQ2 



2 / ArM2 



S_voli 



^1 S_V0l2 



(48) 



/^Qshear q 



N (u\ — U2^^ 



S_V0l2 . 



(49) 



This formalism for the viscosity in the radial direction closely re- 
sembles - except of course for the s hear part - th e c ustomary ID 
viscosit y description given, e.g., bv lDorfil (119981) or lFeuchtinge"J 
([T999ah . 

In the horizontal direction, discretization yields a viscous 
force 



r j 3/A,Mfl 

r [ \ A,.r 

4A? Mfl „ , 
J r 

and a viscous energy dissipation 



Ug 

r 



(50) 



iQSl 



-/^Qesheai- 



Arr 



y I S_voli 



iQffi - -/^QSstiear I -^-y - y I S-VOI2 

4A^ «fl ^ , 



(51) 



(52) 



4.2. Radiative transport 

In the moment description of radiation, the second moment of 
the intensity - which corresponds to the radiation pressure - is 
assumed to be of the following form 



fKr, 



ee 



Km , 



Kn 



J-K,,. 



J-K„. 
2 ^ 



(53) 



with the radial component given by a scalar Eddington factor 

Krr = UaaJ ■ (54) 

The same Eddington factor is taken for both columns 

Krr,\ - feddJ\ Krr,l - feddJl , (55) 

and, for the examples presented in this paper, it has been set to a 
constant value of /edd = 1/ 3 for simplicity. 

The discrete equation of radiative flux in the horizontal direc- 
tion (i.e. for Hg), does not use the full time-dependent equation 
Eq. [3T|but only its stationary part 



V ■ K -H K^pH = 



(56) 



In this way, we did not have to discretize the horizontal compo- 
nent of H ■ Vm, which is, in a similar way to artificial viscos- 
ity, ambiguous in interpretation in the context of the 2C-scheme. 
Considering the simplistic discretization of horizontal exchange 
between the two columns, this stationary, diffusion-like descrip- 
tion remains sufficient. 

Note that the assumption for the radiative pressure in Eq. |53] 
will in general not be consistent with the horizontal radiative 
flux computed from Eq.|56] However, a more consistent descrip- 
tion is not reasonably possible given the limited resolution in 
horizontal direction of the 2C-scheme. Adopting a more elabo- 
rate description would also require solving the detailed 2D ra- 
diative transport to obtain the required Eddington factors (e.g. 
Krr = /edd;-;- J and Kgg = /eddflfl J)- And after all, there is no point 
in improving the radiative transport beyond the level of approx- 
imation of the hydrodynamics part. 

4.3. The discrete equations of radiation liydrodynamics 

Using the discretization scheme presented in Sect. |2] and the re- 

we obtain the discrete version 
3]provides the full discrete set 
of equations of radiation hydrodynamics. These physical equa- 
tions are completed by the equations for the radially averaged 
densities, Eqs.|60]&|6Tl and by the adaptive grid equation. 

As an example of the discrete form of conservative equations 
and to illustrate the notation of the advective contributions, the 
discrete equations of continuity are given here for both columns 



suits from Sect. |4j] and Sect. 143 
of Eqs. |26]-[3l]& Eq. |56] Table 



6[pi S_voli] + 



piS_fluxi 



piS_fluxi 



-Hp'iS_fluxe = 0(57) 
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5[p2S-V0l2] -H p2S_fluX2 -p2S_fluX2 -p2S-fluXfl = .(58) 

L 1/+1 L li 

The notation d[X] indicates a difference of X between the new 
and old time level separated by the time step 6t. Advection rep- 
resented by the terms with SJlux, occurs in the radial direction 
both at the radii r,+i and r, as well as in the horizontal direc- 
tion over the interface between the two columns. Note that the 
horizontal transport terms with SJluxg coiTespond to each other. 

In Table [3] an abbreviated notation was adopted for the ad- 
vective terms by summarizing all three contributions. In this 
form, advective terms, e.g. those from Eq.|57l are written as 



^pi S_fluxi,e s 



pi S_fluxi 



pi S_fluxi 



+pi S_fluxe.(59) 



Spatial differences are denoted as A,- and Ag in the radial and 
horizontal direction, respectively. Averaged quantities - where 
the precise definition depends on the context - are written with 
overhead dashes. 

4.4. The stencil 

The discretization of the system of differential equations at the 
grid point r, also incorporates variables from adjacent grid loca- 
tions. In the present case, dependencies are included up to a dis- 
tance of two grid points. Equations at the grid point / may there- 
fore include variables from / - 2, ; - 1, /, / -H 1, / + 2. Accordingly, 
this ensemble of five grid points is referred to as '5-point stencil' . 

The shape of the stencil is coiTelated closely with the implicit 
solution method because it determines the structure of non-zero 
entries in the Jacobi matrix. In the present implementation, the 
Jacobian is constructed 'ID-style', i.e. all variables from both 
columns (as assembled in Table [1) have only one running index, 
the radial grid point index /. Alternatively, it would also be possi- 
ble to use two running indices as in a 2D code, the second having 
values of only 1 and 2 to differentiate between the two columns. 
This type of indexing would assign fewer variables, only those 
from one column, to each pair of indices, but correspondingly 
also involve a larger 5x2 stencil and a significantly more com- 
plicated algorithm. For 2D grids, this results in a Jacobian (com- 
posed of more numerous but smaller submatrices) that enables 
an increase o f up to 50% in the speed of the matrix inversion 
(IStoklLl2006h . However, in the present (extreme) case, where the 
grid has just two grid points in one direction, the inversion time 
is almost identical to the far simpler ID-like discretization. 

4.5. Averaged density p 

To develop an expression for the momentum in the radial direc- 
tion for the equation of motion, the (scalar) densities must be 
averaged for the same (vector) localization of the velocities (see 
Fig[T]). As a second order advection scheme is used in the radial 
direction, the momentum - and consequently the averaged den- 
sity - is required at 5 successive radius points. Averaging for 5 
successive points is not possible within the 5-point stencil, and 
therefore an additional variable, the radially averaged density p, 
was introduced 



Pi = 



Pi 



(S_volipi|.-i-S_volipi|. J 



V_voli 

i (S_V0l2P2|,. + S_V0l2P2|,._j) 



V_V0l2 



(61) 



4.6. Boundary conditions 

Two successive ghost cells - corresponding to the 5-point dis- 
cretization - constitute the boundary conditions in each column 
at both the inner and outer boundary. 

The inner boundary conditions are stated at a fixed inner ra- 
dius of the computational domain and characterized by constant 
values for p, e, m, J, and H (the same in both columns). The 
value of H entering at the inner boundary corresponds to the lu- 
minosity of the modelled star; m is the mass of the central core. 
The radial velocities at the inner boundary are taken to be zero. 

For time-dependent computations of stellar pulsations - the 
principle task to be solved by the code - the entire envelope of 
the star must be considered. Since nuclear energy generation is 
not implemented in the code, it is impossible, however, to model 
the stellar core. Therefore, the inner boundary is usually placed 
as deep as possible, while remaining clear of the core region 
where nuclear burning might occur. In the case of the Cepheid 
models presented in Sect.|6l the radius of the inner boundary was 
set to be 10% of the photosphere radius. 

The outer boundary conditions are defined at the outermost 
grid point, which moves in a Lagrangian manner, i.e. there is 
no fluid flow over the outer boundary to, or from, the exterior 
space. Accordingly, the radial velocities in the two columns are 
required to be identical at the outer boundary. A common equa- 
tion of motion, formed as the sum by the individual equations 
of motions, determines the gas velocity at the outermost grid 
point - and by means of the Lagrange condition - the veloc- 
ity of the grid point itself. This setup has the advantage that the 
outer boundary of the grid can follow radius variations of the star 
e.g. due to stellar pulsations or structural resettling. Obviously, 
there is no convective flux over the outer boundary. The location 
of the outer boundary in relation to the mass structure is deter- 
mined from the initial model and usually given by a predefined 
ratio (e.g. 1/100) between gas pressure at the outer boundary and 
the photospheric gas pressure. 

For the physical conditions in exterior space, which affect 
the common equation of motion at the outermost grid point, 
% = %■ - ^ ^nd Q = are assumed. These boundary conditions 
are, however, by no means unique and e.g. pext = const, and 
gext = const., or = would also be appropriate. When stel- 
lar pulsations are considered, these outer boundary conditions 
become more influential as they affect the wave reflection and 
dissipation properties. 

The boundary conditions for the radiation field assume free 
radiation at the outer boundary; H is then computed to be // = 
ju 7, where p = | in the case of the Eddington approximation 

/edd - \, and y is a radially averaged value of J. This condition 
is evaluated individually for both columns, so that, in general, 
there will be a different radiative flux from each column. 



4. 7. Temporai centering 

The system of equations of radiation hydrodynamics consists 
of parabolic differential equations. Splitting them into a time 
derivative and spatial terms, they can be written in the form of 



(60) dX 
dt 



^H(X) 



(62) 



These algebraic equations, Eq. |60] & |6T1 are solved implicitly 
together with the system of discrete equations given in Table [3] 



where H{X) is a nonlinear spatial difference operator. The time 
derivative is discretized to be 5[Z]/(5f where 6t is the time step, 
and 5{X] represents a difference in time between the new and 
old time level. To achieve (almost) second order accuracy in 
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Equation of continuity 



Integrated mass (Poisson equation) 
Equation of motion - radial direction 
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5|p,S.vol,] + 


Pi 


S_fluxi| 


L-l 


Pi 


SJluxij 


6[p2 S.V0I2] -1- 1 


P2 


S-fluxzj 


L-l 


P2 


SJluxjj 



A^m = pi S.voli + p2 S_vol2 



(P2M2) V_flux2 + c/247r r A,-(P2) St h — — pi V_vol2 6t /C2P2 Hi V_vol2 dt + Eq2 St = 

r'- c 

Equation of motion — horizontal direction 



Gm _ 
Gm . 



— , um _ An 

(p,Mi) V_fluxi,e + cf\Anr Ay{P{)St + — p; V_voli St ^ipi H\ V_voli St + Fqi <5/ = 



c 

An_ 

c 



(5|peMe H.vol] + (paMfl) H_flux + ^(^2 -^'i) " — (f/i'^iPi + cfiKiPi) HgYl.yoXSt 1- +Fq8& = 

Equation of energy 

Slptei S.voli] -I- ^ (piei) S.fluxi.fl + Pi (c-/i4;rAr(r^iii) + AinifWe) ft - 4;rA-ipi (7] - 5 i)S_voli St - Eanhi - £'adv + ^qi St + Eqoi St = 

S\p2e2 S.V0I2] + (Piei) S.fluX2,fl + Pi (cf2A7lAr(r^U2) - AmifUg^ St - A7IK2 p2 Ul - Si) S_V0l2 St - £anhl - ^adv + £q2 St + Eqoi St = 

S[Ji S_voli] + ^Ji S_fluxi_s + c (c/i47rAr(r^//i) + Ai„,f//a) St+ 



Radiation energy equation 



1 X '^l J \ ~ I^rr 1 

-l-c/i A:,-,-i 47rA,-(rKi)(5f + (7i - 3A:,-,-i) — S.voli St + — -A^^ifUsSt + c ki pi (Ji - 5i)S_voli St = 

r 2 

S\]2 S-V0I2] + 2 -/2 S.fluX2.e + C {cf2ATxAA/H2) - Ai„,f//fl) 5f+ 



-l-c/2 ^rr2 47rAr(r^M2) (5f + (7t - 3A:rr2) S.V0I2 St ■ 

r 



■/2 - '^)7,2 



Aintf Mfl (5r + CK2pl(Jl - S 2) S-VoIt = 



Radiation flux equation — radial direction 

S[Hi V.voli] + V_fluxKe + c/i c47rr2A,(Ji:,,.,) & + f 



3K,-,-i — J\ 

'■ V_voli St + cfi Anr-Hi A^^i) St + c Kipi Hi V.voli St = 



— 1 3K,-,-2 ~~ Jo 1 

H2 V_fluX2,e + Cf2 C Anr Ar(Krr.l) St + C '■ V_V0l2 St + f /2 AnrH2 Ar(U2) St + C A-2P2 H2 V_V0l2 (5/ = 



Radiation flux equation — horizontal direction: stationary limit 

Aiai \Jl~ ^n-.l J 1 - ^rr.l 



I + {cfKip, + cf2K2pl) //eH_VOl = 



time, the spatial terms must be evaluated centered in time, i.e. 
at a point in-between those two time levels in the temporal dif- 
ference. This centering is completed in terms of variables, i.e. 
in the form of //(Zcem.) with Zcem. = 1/2 (x™" + X°^'^), as op- 
posed to centering the operator H^sntiX). This centering of vari- 
ables usuall y provi des a higher temporal accuracy of the scheme 
dPorfi et al.L l2006l) . Based on the centered primary variables, 
successively all other required variables and expressions, such 
as cell volumes, advection fluxes, viscosity terms, opacities, and 
equation of state can be assembled. 



5. Method of solution 

The system of nonlinear, discrete equations is solved time- 
dependently using an implicit Newton-Raphson iteration. The 
implicit solution has the ad vantage of not being affected by the 
CFL time step limit (after ICourant. Friedrichs & LewM, Il928l) 



and also allows the inclusion of elliptical parts into the system 
of physical equations (Poisson and grid equation). The long time 
steps that are possible with the implicit scheme are particularly 
useful for the present problem of convective transport because 
they permit a rapid progression towards the stationary solution. 

Each step in the Newton-Raphson iteration requires the in- 
version of the Jacobi matrix. According to the system of 16 equa- 
tions (Table[T]), the Jacobian is composed of 16 x 16 submatrices, 
which form a pentadiagonal structure of non-zeros reflecting the 
discretization with a 5-point stencil. The inversion of the Jacobi 
matrix uses the customary approach of a Newton-elimination of 
the two lower sub-diagonals, followed by a back substitution 
of the resulting upper triangular matrix. Normalization of the 
Jacobian prior inversion, using the largest term in each discrete 
equation, significantly improves its numerical properties. 

The time step 6t used for advancing the system of physical 
equations is regulated to maintain reasonable iteration numbers 
(usually between 2 and 4) and according to other requirements. 
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e.g. limiting the relative changes in the primary variables per 
time step. In case of divergences, the Newton-Raphson iteration 
is restarted with a reduced time step. 

The crucial point about implicit methods is the computa- 
tion of the derivatives required for the Jacobian. Derivation of 
the discrete physical equations (Table |3]l with respect to the pri- 
mary variables leads to rather elaborate expressions. The implicit 
scheme is also very sensitive to errors and inaccuracies in these 
derivatives. Computer algebra was therefore adopted to allow a 
fast and reliable computation of all required derivatives. These 
computer algebra scripts directly produce FORTRAN code that 
can be plugged-in into a source code. This feature proved to be 
very useful at the development stage because it facilitated nu- 
merous and quick tests of the discretization scheme. 

The computing time for inversion of the Jacobi matrix scales 
with npx ng^ with ng the number of equations (here 16) and np 
the number of grid points (usually np - 500). A current CPU at 
3 GHz achieves about 10 iteration cycles (i.e. time steps) per sec- 
ond for this setup. Unfortunately, the inversion of the Jacobian 
does not parallelize efficiently. Nonetheless, the long time steps 
possible with the implicit solution method ensure that the 2C- 
scheme is much faster than 'classical' explicit 2D or 3D hydro- 
dynamics. 



6. Demonstrating example 



T 




-10 - 



-15 L I I I I I I = 

31 32 33 34 35 36 

Radius [Rg] 

Fig. 6. Details of a Cepheid convection zone: The upper panel 
shows the convective transport in units of the total luminosity 
(solid line), the temperature difference between up- and down- 
drafts (dotted line) and the run of the entropy through the model 
(dashed line, without scale). The convective velocities are given 
in the lower panel: updraft (dashed line), downdraft (dotted line) 
and horizontal (solid line). A positive sign of the horizontal ve- 
locity corresponds to a flow from column 1 to column 2, i.e. from 
updraft to downdraft. The figure focuses only on the outer con- 
vective region, the model actually extends down to about 2>.6Rq. 



According to the intention of applying this scheme in com- 
putations of Cepheid pulsations, a typical Cepheid with Teff = 
5AQQK, L = 10^ Lq, and M = 4.75 Mq (which translates into 
^phot = 36.1 Rq and \ogg - 2) was adopted for testing. This star, 
with a comparatively weak and shallow photospheric convection 
zone, has the advantage that it allows starting from a purely ra- 
diative initial model. For stars with fully convective envelopes, 
this is no longer possible because a purely radiative stratification 
would be too far off and therefore cause a violent collapse of 
the envelope with the onset of convection. The inner boundary 
of the models was placed at 10% of the photospheric radius (i.e. 
~3.6 Rq), although subsequent figures only indicate the outer 
convective region of interest. The models consist of 500 radial 
grid points, the majority of which, due to the adaptive grid, clus- 
ter around the photosphere and in the convective region. 

To model a convection zone with wide up- and more nar- 
rowly confined downdrafts, updrafts are (arbitrarily) assigned 
to column 1 and the corresponding relative cross-section c/i 
is set to a value above 1/2. Accordingly, column 2 covers a 
smaller cross-section and contains the downdraft flows. To en- 
sure that convection finally occures in the intended sense of ro- 
tation, the initial model is perturbed with small radial velocities 
(m < 1 m/sec) using the Schwarzschild convection criterion as a 
guide. 

Starting the time-dependent simulation from that initial 
model, the convective velocity field develops rapidly and grows 
downwards from the photosphere. After a dynamic phase of 
growth that lasts about a thermal timescale of the relevant part of 
the envelope (~10^ seconds), the convective velocities approach 
a stationary solution. The time step then increases quickly and 
the computation is terminated at an age of lO'^ seconds. This 
evolution typically takes around a minute on a 3 GHz CPU and 
requires about 1000 time steps that increase in length during the 
computation from a few seconds at the start up to lO" seconds 
for the stationary solution. 

Figure |6] shows the resulting convection zone using = 
995 1 convective cells on a sphere, donor cell advection for the 
horizontal transport (Eq. [18] A = 0), and a downstream cross- 
section of 20% of the sphere (c/i = 0.8). As a useful guide, 
one can estimate the horizontal scale of photospheric convec- 
tion (in the 2 C-scheme, th i s corr esponds to D, Eq. [TJ to be 
about lOHpQ jFrevtag et al.Lll997h . where Hpo is the character- 
istic photospheric pressure scale height, HpQ - HTsff/g. For the 
present example, a length scale of 20 Hjm was adopted, which 
(evaluated at the photospheric radius) translates into the afore- 
mentioned odd number of cells A^ - 995 1 . This set of param- 
eters, D = IQHpQ, cfi - 0.8, A - Q, serves subsequently as a 
reference for exploring the influence of the individual parame- 
ters. 

The convection zone in Fig.|6]includes the H/He I as well as 
the He II ionization zone. Both are apparent in the entropy pro- 
file given in the upper panel, the former causing the steep photo- 
spheric drop, the latter appearing as moderate gradient between 
34 and 35 Rq. The continued gradient in-between those two ion- 
ization zones (i.e. outwards about 35 Rq) is an effect of the con- 
vective transport and not present in purely radiative models. The 
large temperature difference between up- and downdrafts in the 
outer part reflects a different radial position of the photosphere 
in the two columns. 

The convective velocities (Fig. |6] lower panel) show a com- 
paratively slow updraft motion. In the thin outer regions - around 
the photosphere - the hot material loses energy by radiation 
and changes over to the downdraft column. Due to the narrower 
downdrafts, the downward velocity is accordingly higher and the 
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large momentum in the downdraft motion produces a prominent 
inward overshoot. The mild entropy gradient in that part of the 
envelope also offers only little resistance to the downdrafts. The 
temperature difference is reversed in the overshoot; the down- 
stream flow is now hotter than the 'surroundings', and the con- 
vective flux has a negative sign. Because of the contribution of 
the kinetic energy flux to the convective transport (see Fig. |7]i, 
the temperature difference and the convective flux do not change 
their sign at exactly the same depth. In the overshoot region, the 
material also returns to the updraft column, closing the circulat- 
ing convective motion. 



T 1 1 T 




J I I L 



31 32 33 34 35 36 
Radius [Rq] 

Fig. 7. Contributions to the convective energy flux: transport of 
internal energy by fluid motion (dashed line); kinetic energy flux 
(dash-dotted line); and flux due to work against gas, viscous, 
and radiative pressure (dotted line). The radiative flux makes up 
for the difference between the sum of these species (solid line, 
commonly referred to as 'convective energy flux') and 100%. 



Figure |7] indicates the contributions to the convective flux: 
transport of internal energy, kinetic energy flux, and flux related 
to work against the total pressure (consisting of gas, viscous, and 
radiation pressure). Even though viscous and radiative pressure 
have been included for completeness, the total pressure for this 
type of star is dominated largely by the gas pressure. Radiation 
pressure accounts for up to about 15% of the total pressure, vis- 
cous pressure for much less. Transport of potential energy is not 
evident in Fig. |2l since the contributions from up- and down- 
drafts balance each other in the stationary case. The flux of ki- 
netic energy is entirely inward because of the narrower and more 
rapid downdrafts, which transport more kinetic energy than the 
updrafts. This behavior of the kinetic energy flux is consistent 
with the results from multi-dimensional simulations of convec- 
tion. 

The ability of the code to reproduce the kinetic energy 
flux as well as the extended lower overshoot in qualita- 
tive agreement with multi-dimensional hydrodynamics com- 
putations (R oxburgh & Simmons, 1993; Muthsam et al., 1 19951 : 
ISteffen etal.l i2005) is an indication that the 2C-scheme succeeds 
in describing the essential physics of convective transport. 
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Fig. 8. Convective flux assuming the typical horizontal length 
scale D to equal 2 (dotted line), 10 (dashed line), 20 (solid line), 
and 80 (dash-dotted line) times the characteristic photospheric 
pressure scale height Hj,o. 



6.1. Parameter studies 

The effect of the typical horizontal length scale D - which trans- 
lates into a certain number of convective cells on the sphere 
- on the convective transport is shown in Fig. [8] For values 
of D between IQHpo and SOHpo, the convective flux is only 
slightly affected although the convective flux is somewhat lower 
for large convective cells, especially in the overshoot region. For 
even larger cells, it becomes increasingly more difficult for the 
convective circulation to bridge the growing distance between 
up- and downdrafts, and convection finally ceases. At the other 
extreme, convection also becomes less effective for convective 
cells smaller than IQHpQ. This seems reasonable as many thin 
downdrafts will dissolve rapidly, whereas a smaller number of 
more massive downdrafts can retain their downward momentum 
much longer. This causes the H/He I and He II convection zones 
to separate, as is already apparent in Fig. [8] for the convective 
flux for D - lOi/po (dashed line). Ultimately, there remains only 
a narrow convective region related to H ionization as shown by 
the dotted line for D - 2//po- 

This transition from a large common convective region con- 
taining both the H/He I and He II ionization zone to two de- 
coupled convective shells also happens in a sequence of mod- 
els when changing to 'less-convective' stellar parameters (e.g. 
higher effective temperature). Usually - at least for the inves- 
tigated Cepheid-like stars - the inner He II convection carries 
only marginal flux, although showing convective velocities of 
a several km/sec. These decoupled convection zones found for 
hotter Cepheids are similar to those obtained for A-type stars 
jKupka & Montgomervl 120021; ISteffen et all l2005h . The differ- 
ence in effective temperature of about 2000 K between hot 
Cepheids and cool A-type stars appears to be largely compen- 
sated by the higher surface gravity of the A-type stars. 

In Fig.|9l the effect of cfi is studied by showing convection 
zones with downdrafts taking 25%, 20%, and 15% of the sphere. 
Note that more narrow downdrafts lead, due to correspondingly 
more rapid downdraft motion, to a more pronounced overshoot 
despite a reduced overall convective efficiency. Concerning the 
efficiency of convection, a 50/50 ratio of up- and downstream 
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Fig. 9. Influence of the ratio of cross-section between up- and 
downdraft. The plot gives the convective fluxes of convection 
zones where the downdrafts take 25% (dashed line), 20% (solid 
line), and 15% (dotted line) of the sphere (i.e. cf\ is 0.75, 0.8, 
and 0.85). 



Increasing A from to 0.5 changes the horizontal advection 
from donor cell to centered values, which successively reduces 
the dissipation of radial momentum due to horizontal exchange. 
Consequently, the downdraft moment is retained longer when 
the convective circulation makes its turnaround in the lower 
overshoot region. The effect of increasing A is therefore basi- 
cally a deeper overshoot as well as a higher overall convective 
efficiency because of the reduced dissipation. 

Summarizing the discussion of the parameter influence and 
considering qualities such as convective flux, depth of the con- 
vective region, and amount of overshoot, it appears that, al- 
though there is some variability in the results, the basic type of 
solution is quite robust. It is possible to suppress convection by 
choosing extreme parameters, some combinations of parameters 
may also cause numerical problems; none of the test computa- 
tions, however, produced a convective region qualitatively dif- 
ferent from those shown in Figs. l8]-fT0l 

Even though the adopted parameters are up to now little 
more than an educated guess and still require verification by 
comparison with observations or more elaborate numerical sim- 
ulations, the 'reasonable parameter range' suggested here is 
probably quite reliable. 

6.2. Accuracy of the discretization 



cross-section would obviously be the optimum, but that is prob- 
ably not a realistic scenario for photospheric convection in real 
stars. In contrast to the horizontal length scale D, for which the 
hydrostatic pressure scale height //po provides good indications 
of a reasonable parameter range, the proper value of cf\ is more 
difficult to estimate and requires further investigation. The gran- 
ulation pattern of the Sun as well as multi -dimensional hydro- 
dynamics computa tions of other stars (e.g. lFrevtag et al.Lll996l; 
ISteff^en et al.Ll2005l) clearly suggest rather narrow downdrafts. 
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Fig. 11. Accuracy of flux conservation throughout the convec- 
tion zone presented in Fig.|6l For stationary solutions, deviations 
from a constant energy flux correspond to errors in the total en- 
ergy conservation. Inwards of the plotted region, the flux is en- 
tirely transported by radiation and hence no errors occur. The 
same is true for the outermost part of the model above the con- 
vective region. 



Fig. 10. Effect of the parameter A for the horizontal advection 
(see Eq. [TST l. The figure shows the relative convective flux for 
A equalling (i.e. donor cell, solid line), 0.05 (dashed line), 0.1 
(dotted line), and 0.15 (dash-dotted fine). 

Figure [10] illustrates the effect of the horizontal advection 
scheme quantified by the parameter A as described in Sect. 12.41 
As already argued there, one should keep A well below 0.5. 



According to the nature of convection, a considerable part of 
the luminosity is converted from radiation to internal and kinetic 
energy, transported upwards through the convection zone, where 
radiation once again takes over. This conversion of energy causes 
errors in the total energy balance, which are eventually evident in 
the total energy flux (i.e. luminosity) for stationary solutions. In 
the present type of discretization, the total energy is not treated 
conservatively but is composed of several species (internal, ki- 
netic, potential, and radiative energy); the total energy conserva- 
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tion is therefore a good measure of the discretization accuracy 
and a possible way of testing physical soundness. 

The deviations from constant total energy flux are shown in 
Fig.[TT]for the standard-parameter convection zone. Without the 
correction term accounting for the momentum dissipated by hor- 
izontal advection given by Eq. |24l the discrepancies would be- 
come larger than 30%. 

7. Conclusions 

The scheme proposed in this paper has a number of advantages: 

- It is a non-local description of convection and therefore pro- 
vides a consistent computation of the depth of the convective 
region including convective overshoot. 

- The convective flux is computed directly from hydrodynam- 
ics and not from a heuristic, parametrized model. 

- The two-column convection has stationary solutions and in 
principle allows arbitrarily large time steps. This implies that 
it is suitable for application to problems involving long time 
series, such as stellar pulsations or stellar evolution. 

- The 2C-model is much faster than multi-dimensional hydro- 
dynamics computations; stationary solutions can be obtained 
within minutes. 

- Radiative transport is an intrinsic part of the scheme, i.e. no 
hydrodynamical model with plugged-in radiation effects. 

- The main parameters of the scheme have a straightforward 
geometrical meaning that also provides indications of rea- 
sonable values for these parameters. 

However, there are also shortcomings to be considered: 

- The 2C-scheme is basically still a parameter-dependent 
model. These parameters require proper adjustment. 

- Horizontal advection and radiative transport are poorly rep- 
resented because of the very coarse 'two-cell' spatial resolu- 
tion in the horizontal direction. 

- The 2C-model uses a simplistic description of the full spec- 
trum of vertical and horizontal convective motion, which ig- 
nores turbulence effects and limits the investigation of more 
subtle features of convection. 

In its present form, the two-column scheme provides a sim- 
ple, yet physically sound and consistent, non-local, radiation- 
hydrodynamics description of the convective circulation. The 
model still contains free parameters, but their geometrical in- 
terpretation provides at least reasonable indications of proper 
values, and they do not change the results by magnitudes. For 
applications in which more detail and higher certainty is re- 
quired, more elaborate methods, such as multi-dimensional hy- 
drodynamics or turbulence models, remain the most appropriate 
alternative. 
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